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AN UNCONDITIONALLY STABLE DISCONTINUOUS GALERKIN 
METHOD FOR THE ELASTIC HELMHOLTZ EQUATIONS WITH 

LARGE FREQUENCY 

XIAOBING FENG* AND CODY LORTONt 

Abstract. In this paper we propose and analyze an interior penalty discontinuous Galerkin 
(IP-DG) method using piecewise linear polynomials for the elastic Helmholtz equations with the 
first order absorbing boundary condition. It is proved that the sesquilinear form for the problem 
satisfies a generalized weak coercivity property, which immediately infers a stability estimate for the 
solution of the differential problem in all frequency regimes. It is also proved that the proposed IP- 
DG method is unconditionally stable with respect to both frequency uj and mesh size h. Sub-optimal 
order (with respect to h ) error estimates in the broken H 1 - norm and in the L 2 -norm are obtained 
in all mesh regimes. These estimate improve to optimal order when the mesh size h is restricted to 
the pre-asymptotic regime (i.e., = 0(1) for some 1 < f3 < 2). The novelties of the proposed IP- 

DG method include: first, the method penalizes not only the jumps of the function values across the 
element edges but also the jumps of the normal derivatives; second, the penalty parameters are taken 
as complex numbers with positive imaginary parts. In order to establish the desired unconditional 
stability estimate for the numerical solution, the main idea is to exploit a (simple) property of linear 
functions to overcome the main difficulty caused by non-Hermitian nature and strong indefiniteness 
of the Helmholtz-type problem. The error estimate is then derived using a nonstandard technique 
adapted from [9]. Numerical experiments are also presented to validate the theoretical results and 
to numerically examine the pollution effect (with respect to lj) in the error bounds. 

Key words. Elastic Helmholtz equations, Korn’s inequality, unconditional stability, discontin¬ 
uous Galerkin methods, generalized weak coercivity, error estimates. 
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1. Introduction. Wave phenomena are pervasive throughout many scientific 
fields, wave computation has been one of the central topics in computational sciences 
and has also been at the forefront of scientific computing. It becomes more and 
more important as the boundary of wave-related application problems keeps pushing 
outward. Solving these application problems largely hinges on computing the solutions 
of the governing wave equation(s). Among them those involving high frequency waves 
are often most difficult to analyze and solve numerically. This is due to the fact that a 
high frequency (or large wave number), which means a very short wave length, results 
in a strongly indefinite PDE problem and a highly oscillatory solution. Consequently, 
very fine meshes must be used to resolve the wave, which in turn results in extremely 
large non-Hermitian and strongly indefinite algebraic systems to solve, a daunting 
task especially in high dimensions. 

In this paper we shall focus our attention on one type of wave phenomena, that is 
wave propagation through some isotropic elastic media. This type is often encountered 
in applications from materials science, medical science, petroleum engineering, just to 
name a few. Specifically, we shall consider the following elastic Helmholtz problem: 

(1.1) — uj 2 pu — div (cr(u)) = f in H, 

(1.2) \uiAu + <r(u)n = g on I := dfl, 
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where 12 C M. d ,d = 2,3, which represents an isotropic elastic media, is an open 
and bounded domain, n denotes the unit outward normal vector to 912, and i = 
\J — 1 denotes the imaginary unit, u : 12 —► C d denotes the (Fourier-transformed) 
displacement vector of the elastic media 12. ui and p are both positive constants 
for which u denotes the frequency of the sought-after elastic wave traveling through 
12 with speed p. A is a d x d real symmetric positive definite matrix. <j(u) is the 
stress-tensor in 12 defined by 


<t(u) := 2/ue(u) + Adivu/, 


e(u) :=i(Vu + Vu T ) 


where p, > 0 and A > 0 are the Lame constants for the elastic media 12. 


Equations (1.1) and (1.2) can arise either from the frequency domain treatment 


of linear elastic wave equations (cf. [2]) or from seeking a time-harmonic solution to 
the following linear elastic wave equations: 


(1.3) 

(1.4) 

(1.5) 


pV tt - div (cr(U)) = F 
AU + cr(U)n = G 

u = U t = 0 


in 12 x (0, oo), 
on r x (0, oo), 
in 12 x {t = 0}. 


When G = 0 the boundary condition (1.4) is known as the first-order absorbing 


(radiation) boundary condition [5j. This boundary condition is an artificial boundary 
condition which has the property of completely absorbing incoming plane waves that 
are perpendicular to the boundary 912. This boundary condition is imposed on the 
boundary 912 of the truncated domain 12 of some unbounded domain with the intent 
of making the problem computationally feasible. 

As is the case for the scalar Helmholtz equation, which governs time-harmonic 
acoustic waves, the frequency to also plays a key role in the analysis and implementa¬ 
tion of any numerical method for the elastic Helmholtz equations. It is a well-known 
fact that in order to resolve the wave numerically one must use some minimum num¬ 
ber of grid points in each wave length t = 2n/io in every coordinate direction. This 
requires a minimum mesh constraint toh = 0(1), where h denotes the mesh size. In 
fact, the widely held “rule of thumb” is to use 6 — 10 grid points per wave length. In 
m, Babuska et al proved the necessity of this “rule of thumb” in the 1-dimensional 
case of the scalar Helmholtz equation. In p2] it was also shown that the H 1 error 
bound for the finite element solution contains a pollution term that contributes to 
the loss of stability for the finite element method in the case of a large frequency. 
This pollution term also causes the error to increase as u increases under the mesh 
constraint uh = 0(1). This forces one to adopt a more stringent mesh condition to 
guarantee an accurate numerical solution. 

The loss of stability of the standard finite element method applied to Helmholtz- 
type problems remains as an open issue to be resolved. In 1110 a stringent mesh 
condition of co 2 h = 0(1) (called the asymptotic mesh regime) was required to obtain 
optimal and quasi-optimal error estimates for the scalar Helmholtz equation. In |j a 
similar mesh constraint of the form ui^h = 0(1) was used to obtain error estimates 
for the elastic Helmholtz equations, where /3 > 2 is a constant to be discussed later in 
this paper. Requiring such a stringent mesh constraint makes the use of a practical 
coarse approximation space impossible in the case that uj is large. This is a major 
hurdle that must be overcome if one wishes to use multi-level algebraic solvers such 
as the multigrid method or the multi-level domain decomposition method. 
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The primary goal of this paper is to develop and analyze an interior penalty 
discontinuous Galerkin approximation method which is unconditionally stable with 
respect to both u and h for the above elastic Helmholtz problem, which is different 
from the goal of the earlier DG work M- Specifically, we seek a numerical method in 
which a priori solution (or stability ) estimates can be obtained for all w, h > 0. Our 
work in this paper can be regarded as the elastic counterpart of those reported in mm 
m, where the scalar Helmholtz equation and the time-harmonic Maxwell’s equations 
were the focuses of study. We like to note that all three Helmholtz-type problems 
share the above mentioned two main difficulties, namely, strong indefiniteness and 
non-Hermitian character. However, due to the fact that to their leading operators are 
very different, in particular, they have quite different kernel spaces, the formulation 
and analysis of the desired IP-DG methods are significantly different, although they 
are conceptually similar. 

The remainder of the paper is organized as follows. In Section [2] we establish 
a new generalized weak coercivity property for the elastic operator —div(<r(-)), and 
derive a stability estimate for the PDE solution as a corollary of the generalized 
weak coercivity property. In Section [3] we introduce some notation and present the 
formulation of our IP-DG method studied in the rest of the paper. In Section [4] we 
derive a stability (or a priori solution) estimate for the numerical solution in all mesh 
regimes including the pre-asymptotic regime (i.e. = 0(1) for some 1 < /3 < 2). 

In Section [5] we employ the non-standard error estimate technique, which uses an 
elliptic projection along with the (discrete) stability estimate, to establish optimal 
order error estimates in both broken H 1 -norm and the L 2 -norm. Section [h] is devoted 
to numerical experiments that validate our theoretical results and demonstrate the 
proposed method’s advantage over the standard finite element method. 


2. Generalized inf-sup condition and stability estimates for PDE solu¬ 
tions. The goal of this section is to prove that the sesquilinear form a(-, •) defined 
in the weak formulation of ( flTlt (l~2f satisfies a generalized weak coercivity prop¬ 
erty. Similar results were also obtained for the scalar Helmholtz equation and the 
time-harmonic Maxwell’s equations mm ■ As a direct consequence of the weak 
coercivity property proved in this section, we will obtain stability estimates for the 
solution u of (1.1)-(1.2). 

Throughout the rest of this paper we adopt the standard L 2 -space norm and inner 
product notation. In particular, for Q C 0 and £ C T, let and (v)e denote 

the L 2 inner product on the complex-valued inner product spaces L 2 (Q) = (L 2 (Q)) d 
and L 2 (£) = (L 2 (£)) d , respectively. 

We will study the weak formulation of (1.1 )—( 1.2 ) defined as seeking u £ H 1 (I2) 
such that 


(2.1) o(u,v) = (f,v) n + (g,v) r Vv £ H 1 (f2), 
where a(-, •) is the sesquilinear form defined on H 1 (H) x H 1 (H) as 

(2.2) a(u,v) := A(div u, div v) n + 2/x(e(u), e(v)) n - uj 2 p(u, v) Q + iw(Au, v) r . 

To show that a(-,-) satisfies a generalized weak coercivity property we will need 
to assume that the domain fl is “nice”. In particular, we assume that Q is a strictly 
star-shaped domain, i.e. there exists Xq £ H and Cq > 0 such that 


(2.3) 


(x — xq) • n > cq Va; £ T. 
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We will need to make use of Korn’s second inequality to achieve our result. This 
inequality is given in the following lemma. 

Lemma 2.1. There exists a positive constant K such that for any v £ H 1 (fl) the 
following inequality holds 


(2.4) 


I hRo) L 


< K 


I £ ( v )IIl 2 (o) + IM| L 2 (n) 


Lastly, we will need to assume that the solution u to (2.1) satisfies a Korn-type 
inequality on the boundary L. This boundary Korn inequality was conjectured in 
0 and it was necessary for the authors to obtain stability estimates for the Elas¬ 
tic Helmholtz equation that are optimal in the frequency w. Here it is stated as a 
conjecture. 

Conjecture 2.2. There exists a positive constant K such that for any u £ 
H 2 (H), the following Korn-type inequality holds 


With this boundary Korn inequality in mind, we define the following function 
spaces. 


V := {v £ H 1 (H); |Vu| r £ L 2 (T)}, 

2.5|) holds with constant I\ 


Vft- ; = {v € V 


}■ 


We also note that since A in (1.2) is SPD, there exist positive constants ca and 
Ca such that 


( 2 . 6 ) 


CA||v||| 2(r) < (dv,v) r < CU||v||| 2(r) Vv £ L 2 (T). 


The following theorem states the generalized weak coercivity property for the 
sesquilinear form a(-, •). 

Theorem 2.3. Let LI be a strictly star-shaped domain with diam{LV) = R and let 
K be a positive constant. Then for any u £ V^> the following inequality holds 


(2.7) 

where 


I Ima(u, v)| 

sup —n—n-1- SU P . 

vgv ll v lls veHqn) lll v ll|L 2 (n) 7 


Rea(u, v)| 1 

> —Hulls, 


7 := 


4 R 2 K (1 + ^1+ 4 R 2 K +(d- l) 2 


2/j 


M := 


^i?w 2 pco/i + R z loCaK + c^p 2 
I E '■= w2 p|l v lll 2 (n) + A||divv||| 2(n) + 2/j||e(v)|| L 2 (n) 


CACo/i 

,2 _ 


+ c 0 A||divv||| 2(r) + c 0 Ml| £ (v)||| 2(r) , 
\l 2 (n) := w 2 p|l v lli 2 (n) + co^lMlisqq. 
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Proof. In this proof we assume that u£ V^fl H 2 (fl), noting that once we obtain 
the result on this more restrictive space that we can obtain the desired result for 
u € ~V K by the standard density argument. 

By letting v = u in (2.2) and taking the real and imaginary parts separately we 
get 

( 2 . 8 ) 

(2.9) 


Rea(u,u) = A||divu|| 2 2(n) + 2^||e(u)|| 2 2(n) - w 2 p||u|| 2 2(n) , 
Ima(u, u) = w(Au, u) r . 


For the remainder of this proof we let v = (Vu)a, where a := x — x 0 and x 0 £ Q 
comes from the star-shape property on the domain given by (2.31. With this choice of 
test function v, we can use the following integral identities for the elastic Helmholtz 
operator [Sj- 


( 2 . 10 ) 


A(a ■ it |divu| 2 ) r + 2p(a • n, |e(u)| 2 ) r 

= (d~ 2)(A||divu|| 2 2(n) + 2/x||e(u)|| 2 2(n) ) 

2Re ^A(divu,divv) n + 2/x(e(u), e(v)) n ^ , 


+ : 


and 

( 2 . 11 ) 


^ll u lll 2 (n) = (a-n, |u| 2 ) r - 2 Re(u, v) n . 


With v = (Vu)a in ( |2.2[ ), we apply the above identities to get 

(u)| 2 ) r 


2Rea(u,v) = A(a • n, |divu| 2 ) r + 2/i( 


a ■ n, \e 


- {d~ 2)(A||divu||| 2(n) +2/i||e(u)||| 2(n) 


+ du 2 p\\u\\ 2 L2{n) - u 2 p(a • n, |u| 2 ) r 
— 2wIm(Au, v) r . 

Rearranging the terms gives 

(2.12) duj 2 p\\u\\ 2 L2{Q) - {d - 2) (A||div u|| 2 2(n) + 2/r||e(u)|| 2 2(n) 

= -A(a • n, |divu| 2 ) r - 2p(a ■ n, |e(u)| 2 ) r + w 2 p( 
+ 2wIm(Au, v) r + 2Rea(u, v). 


a ■ n. u 


2 )x 


Adding (2.12) and (d — 1) times ( |2.8[ ) yields 

u 2 p\M 2 L 2 [Q) + A||divu|| 2 2(n) + 2^||e(u)|| 2 2(n) 

= -\(a • n. |divu| 2 ) r - 2 p(a • n, |er(u)| 2 ) r + uj 2 p(a ■ it |u| 2 ) r 
+ 2uj Im ( Au, v) r + Re a(u, (d — l)u + 2v). 

(6|), and Young’s inequality yields 
2 


Now applying (2.3) 
.2 


w “Pll u IL 2 (n) + A||divu|| L2(n) + 2^|| £ (u)|| i2(n) 

< -c 0 (A||divu||| 2(r) +2/z||£(u)||| 2(r) ) + Rw 2 p||u||| 2(r) 
R 2 uj 2 C a . 


+ 


u llz, 2 (r) + <S||Vu||£ 2(r) + Rea(u, (d — l)u + 2v). 


<5 
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Noting that u £ V^ 


ana cnoosmg 


w 2 p||u|| 

2 

L 2 (f2) 


1 

< 

— 1 


OJCA 


(Roj 

< 



- A||divu||| 2(a) + 2/x||e(u)||£ 2(n) + c 0 (A||divu||| 2(r) + 2/z||e(u)||£ 2(r) ) 
„ o R 2 uj 2 CaK \ „ 2 „ , _ 


9 R 2 oj 2 CaK \ ^ . 

■ 2 p 4-1- Co/iJwc^Hull^p) + Re o(u, (d - l)u + 2v) 


Cod 

• 2 -+ R 

UJCACoP 

M | Ima(u, u 


= | Rea(u, v 
where v := (d — l)u + 2v. Thus 
(2.13) 


cop ' 

uCaK + c Q p \ /^ u \ _|_ Rea(u, (d — l)u + 2 v) 

AC 0 p ) 


u|||; < | Rea(i 


u, v) | + M | Im a(u, u) |. 


By the definitions of || ■ ||^ and ||| • ||| i 2 (n) and using (2.4) we obtain 

\\M\l H n)=^PMl H n) + copM 2 r 

< 4R 2 w 2 p||Vu|| 2 2(n) + 4R 2 c 0 p\\\ 

b (d - l) 2 w 2 p||u|| L 2 (n) + {d- l) 2 c 0 A<||u|| 2 2 (r) 

/ 9 \ -> 


||Vu || 2 


< 


4R 2 K ^1 + j + 4R 2 K + {d- l ) 2 


IM||. 


Thus 
(2.14) 

It follows from ( 2.13| ) and (2.14) that 

|Rea(u, v)| > |Ima(u, u)| |Rea(u, v)| 


^IIU 2 (n) < 7l|u|U. 


I Ima(u, v)| 

SU P-n —u - h SU P 

vGV ||V||_E vGH 1 (0) 


l v IIU 2 (n) II u I|b |||v|||z, 2 (n) 

M | Ima(u, u)| + | Rea(u, v)| 

7IMI-E 

1 


> 


> -l! u IU- 

7 


The proof is complete. □ 

The generalized weak coercivity property of a(-, •) given in Theorem 2.3 immedi¬ 
ately infers the following stability estimate for the solution to ( 2 . 1 ). 

Theorem 2.4. Let Q be a strictly star-shaped domain with diam(Q) = R. Fur¬ 
thermore, suppose there exists some positive constant K such that u £ V ^ solves 
© for f £ L 2 (fl) and g £ L 2 (T). Then the following stability result holds: 

) + II §|| L 2 (r)^ ; 


(2.15) 


IMU<2 7 (4- + —) (ll f IU 2 (0)- 

\u 2 p c 0 pj V 


where || ■ Us and 7 are the same as those defined in Theorem 2.3. 
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Proof. By the Cauchy-Schwarz inequality, for any v £ H 1 (fl) the following result 
holds 




| ) 1 \01 ’ \1 | 

< l|f|U 2 (n)||v||i2 (n ) + ||g||i,2( F ) 

^ (ll f llz, 2 (n) + llslli 2 (r)) (ll v llz, 2 (n) + ll v lli 2 (r) y 

“)lll v IIU 2 (fi)(ll f IU 2 (fi) + ||g||L 2 (r)^ 

-) il v ll-E(ll f l|L 2 (n) + ||g||L 2 (r)), 




UJ 2 P C 0 /i 


— I 2 ' 

Vo rp CQfi; 

Combining this inequality with the result of Theorem |2.3| yields 

|Imo(u,v)| |Rea(u, v)| 

J — ' + sup J v n 


-Hulls < sup -rr— Ji - 1- SUP -T— m - 

7 vev ll v lls veH!(Q) lll v lllL 2 (n) 

- 2 (i + ^)( l|f||i2 ( n ) + l|g||i2 ^)' 


The proof is complete. □ 

Remark 2.5. We remark that both estimates (2.7 1 and (2.15) still hold without 
Conjecture \2.2\ but the dependence of the upper bounds onuj would be slightly stronger 
as shown in m- 


3. Formulation of the IP-DG scheme. In this section we shall present our 
IP-DG scheme for the elastic Helmholtz equations. Conceptually, the formulation 
will mimic those for the scalar Helmholtz problem and the Maxwell problem given in 
nun]. To do so we first need to introduce some notation. 

Let Th be a shape regular triangulation of the domain fl such that for each 
triangle/tetrahedron K £ Th, h > hx := diam(A'). Also for each edge/face e of 
triangle K define h e := diam(e). Next we define the following notation for the set of 
edges of Th' 


(3.1) := set of all interior edges/faces of Th, 

(3.2) := set of all boundary edges/faces of Th on T. 


For each e £ there exists A' e , K' e £ Th such that e = dK e n dK' e and thus we 
define the jump and average operators on e as follows: 



v| k ' , if the global labeling number of K e is greater than that of K' e 
v| k c , if the global labeling number of K' e is greater than that of K e 


and 

Mle := \( v \ k c +v|k')- 

Note that for e £ £}^, we use the convention [v]| e = {v}| e := v| e . Also for e £ £[ we 
shall need to define the outward normal vector n e to e = dK e D <9 AT'. Let n e be the 
unit outward normal vector to K e on e if K e has a greater labeling number than A'', 
and the unit outward normal vector to AT' if the opposite is true. 
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Next, we introduce the energy space E and the sesquilinear form ah{-, •) on E x E 
as follows: 

(3.3) E:= JJ H 2 (A'):= ]J {H\K)) d , 

KGT h KGT h 


(3.4) a^(u, v) := b h ( u, v) + i(j 0 (u, v) + Ji(u, v)) 

where 

bh( u,v):= Y (A(divu,divv)^ + 2/z(e(u), £(v)) K ) 
K&T h 


vY ([ u ]>M v H}) e 


eesl 


Jo(u,v) := Y X^([ U M V ])> 


ee£l 


J i(u, v) := Y 'n,eh e ([a(u)n e ], [<j(v)n e ]). 


ee£l 


V u, v G E, 


Y (M u ) n e}, t V ])e 


Here r/ is an /i-independent “symmetrization” parameter and 70 ,e; and 7i, e are non¬ 
negative “penalty” parameters that will be specified later. For the rest of this paper 
we shall restrict ourselves to the case rj = — 1 (i.e. the symmetric case) but we note 
that other values of 77 are possible. 

Remark 3.1. (a) Both iJo and iJ\ terms are called interior penalty terms. 

To the best our knowledge, penalizing the jumps of the normal stress cr(u)n e across 
the element interfaces seems have not been used before, and this term will play an 
important role in our stability and error analysis to be given in the next two sections. 

(b) Another new feature of the sesquilinear form ah is that the penalty constants 
in iJo and iJ± are complex numbers and must have positive imaginary parts. This 
feature will be vital for constructing our unconditionally stable IP-DG scheme. 

(c) It is easy to see that is a consistent discretization of the elastic operator 

—div (cr(-)). 

We shall also define some appropriate semi-norms/norms on the energy space E. 

1 

Y A H divv lli 2 (K)+ 2 Ml|e(v)|||2 (if) ) , 

IMIl.h : = (l v ll,ft+ Y (■^IlMlI^e) + Ti.e^ell [cr(v)n e ] ||^, 2(e) ) 

IIMIIl ,h ■■= ( l|v||i, h + Y ~ ||M V K} 

The weak formulation for & ( |1.2[ ) based on the sesquilinear form ah(-,-) is then 
defined as seeking u € EnH( oc (Sl) nH^H) such that for all v G EnH^ oc (H)nH 1 (H) 

(3.5) a h ( u,v) - w 2 p(u,v) n + iw(Hu,v) r = (f,v) n + <g,v) r . 
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Our discontinuous Galerkin finite element space is taken as the standard piecewise 
linear polynomial space 

V/i := 1] p i(*0 : = II 

kgTh Ker h 

where I\ (K) is the space of linear polynomials with complex coefficients on the tri¬ 
angle/tetrahedron K £ Th- Then our IP-DG method is defined as seeking 
such that 

(3.6) A h (u h ,v h ) = (f,Vfc) n + (g,v h ) r \/v h £ V hl 
where 

(3.7) A h ( u h ,v h ) := a h (u h ,v h ) - cv 2 p(u h ,v h ) n + icv(Au h ,v h ) r . 

In the following sections we shall establish stability and error estimates for the 
IP-DG method defined above. To prove these estimates we shall require that there 
exist h, 7 q, and 71 > 0 such that 


70,e S 


7i,e < 71 < Cji <e Ve £ £, 


h ? 


for some positive constant C independent of h. It will also be convenient to introduce 
the notation £ := 1 + 7(7 1 - It is easy to show the following norm equivalence on V^: 


(3.8) 


| v h||l,7i < 


l,h < 11V^ 11! ^ 


Vv/,. £ V h . 


4. Stability estimates. The goal of this section is to establish unconditional 
stability estimates for the solution of our IP-DG method, which can be regarded as 
the discrete counterpart of the following frequency-explicit estimate for the solution 


u £ H 2 (fl) of the problem (1.1) (1.2), which was established in [5]. 

Theorem 4.1. Suppose that fl is a convex polygonal domain or a smooth domain 


and u £ H 2 (H) solves (1.1 1 -(1.2). Then u satisfies the following estimate: 


(4-1) IMI-tf 2 (n) < C(w“ + (||f||£,a(ri) + |g||i 2 (r))i 

where a = 1 if there exists K > 0 such that u £ as defined in Section [i| and 
a = 2 otherwise. 

Remark 4.2. The desired Korn-type inequality on the boundary T of a general 
domain tt required for u £ V^> is yet be proved and is posed as a conjecture in | 2 J/, 
although there are strong evidences to indicate that this is the case. Moreover, we 


believe that the optimal dependence on u> (i.e., when a = 1) should hold in (4.1) as 
that is is the case for the scalar Helmholtz equation and the time-Harmonic Maxwell’s 
equations. For the rest of this paper we shall keep this parameter a. 

In f 1 Sf . an a priori error estimate is obtained in the case Lo 1+a h = 0(1) (called 
the asymptotic mesh regime), with stability being a consequence. This analysis is based 
on the so-called Schatz argument m- Similar results in the asymptotic mesh regime 
have been proved for the standard finite element method applied to the scalar and 
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elastic Helmholtz equations (c.f. HGua w- In this paper we are more interested in 
deriving solution estimates in the pre-asymptotic mesh regime, i.e. ui^h = 0 ( 1 ) for 
some 1 < /? < 1 + a. This is the most interesting case because similar results have not 
been proved for the standard finite element method for the elastic Helmholtz problem 
and the analysis requires a non-standard argument. 

We begin by first proving some discrete weak coercivity properties for the sesquilin- 
ear form Ah(-,-)- The desired stability estimates can then be obtained as a conse¬ 
quence of these discrete weak coercivity properties. 

Let us start with a technical lemma which will be instrumental in establishing the 
discrete coercivity of the sesquilinear form We note that this lemma makes 

use of the fact that the space V/, is made up of piecewise linear polynomials, and thus 
a different strategy must be sought in the case that is a higher order DG finite 
element space. 

Lemma 4.3. For any 0 < 5 < 1 and Vh £ there holds 


(4.2) \v h \l h < &u 2 p|K||2 2(n) 

C s 


° S u\\v h \\ 2 L 2 (r) + —J 0 (v h ,v h ) 
7o 


u) 2 ph 2 71 


wh 
Jl(Vfc, Vh) 


where C is a positive constant independent of ui, h, 70 , and 71 . 

Proof. We note that vj,| k £ P 1 (I\) for all K £ Th and thus divcr(v^)|^ = 0. 
For any w h,v h £ and K £ Th, integrating by parts yields 


0 = (div cr(v/,), wj) K 

= -(a(-v h )n K ,w h ) gK + A(divv/ 1 ,divwh) if + 2p,(e(v h ), e(w h )) K . 
Here we have used a differential identity from Lemma 3 of j5], namely 


(4.3) cr(vh) : VWh = Adivv^div w/j + 2p,£{v h ) : e(w h ). 


Thus 


A(divvh, div + 2p(e(v h ), e{w h )) K = (a(v h )n K ,w h ) gK . 

Setting Wh = and summing over all I\ £ Th gives 

\vh\l, h = Y ( a (v h )n Kl v h ) gK 

KeT h 

= Y ((Mv/OM, [vh]) e + <[<7(v h )n e ],{v h }> e ) + Y (o-(vh)n e ,Vh) e 
ee£i ee£Z 

< Y (ll{ Cr (w)n e }|| L 2 (e) ||[Vh]||L2(e) + ||[cr(Vh)n e ]|| L 2 (e) ||{v /l }|| L 2 (e) ) 

£I h 

+ Y ll Cr K)n e || L 2 (e) ||Vh|| L 2 (e)- 
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Applying the trace and inverse inequalities on the right-hand side yields 
Iv„|?, h <C7 £ h~ 2 11 ^ ) 11L 2 (K e ) 11 11L 2 (e) 

+ c 51 ^ 5 IIK]IU 2 (e) (lkK)|| L 2 (ife) + ||<r(v h )|| 1/ 2 (K -, ) ) 

+ C Y1 h e^ \\W{ v h)n e ]\\ L 2 ie) (|K|| L 2 (ife) + llv^ll^^,)) . 


Finally, it follows from the discrete Cauchy-Schwarz inequality along with Young’s 
inequality that 


\vh\l,h < 2 ( 5Z ll cr ( v ^)IU 2 (N) j ||vh||L2(an) 

\KeTh ) 

+ ^7o 2 f 5^ ll[ v ^]lll 2 (e)) ( 5Z ll°'( V ^)ll|2(X) 


e6 £l 


KeT h 

I 2 


C7l *( 5Z 7 i,e/ie||[o-(v/ l )n e ]||i 2(e) )'||v /l || I , 2 (n) 
e£S! 


C{ X + 2p) 


C( A + 2 /i) 


< ^|v h |i )h +-^- w IKIIz7 (9n) + 2 l v *-li,^ +- §l0 J o(v fe ,v h ) 


C 


+ 5(1 - *)« pKH^ + v,). 


Thus (4.2 1 holds. □ 

We are now ready to state the following main theorem of this section. 

Theorem 4.4. Let Ah{-,-) be the sesquilinear form defined in (3.7). Then there 
hold for all e V/,. 


(4.4) I Ah(vh, v/,)| > c(% + + UJ 2 j 1 2 r y 1 ) (l v ^li,/i + w2 pll v ^lli 2 (n)) > 

(4.5) \A h (v h , v, t )| > Jo(v h ,v h ) + Ji{v h ,v h ) +w(4v h ,v h ) r , 

where £ = 1 + and C is a positive constant independent of to, h, 70 , and 71 . 
Proof. Let £ Vfc. Taking the real and imaginary parts of Ah{wh,~Vh) yields 


(4.6) ReA h (v h ,v h ) = \v h \^ h - uj 2 p\\v h \\ 2 L2{n) - 2Re ^ <{cr(v^)n e }, [v h ]) e 


ee£l 


(4.7) Im A h (v h ,v h ) = J 0 (y h ^h) + Ji(v h ,v h ) + u(Av hl v h ) r . 

Then (4.5) follows directly from ( |4.7| ). 

To verify (4.4), we need to bound the term YleeE 1 ({°'( v /») n e}i [ v *.]) e , which in¬ 
volves using the trace and inverse inequality and was already carried out previously 
in the proof of Lemma |4.3[ Thus 


Red/,(v/,,v/,) < \v h \l h - u 2 p\\v h \\ 2 L2{n) +(770 s ( ^l|[vh]||| 2 (fl) ) 2 \v h \i,h 


ee£l 


< 



“ 2 p\\vh\\ 2 L2(n) + 


—Jo(v hl v h ). 

7o 
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Combining the above inequality with (4.21 and using <5 = \ we get 


C 


xM i,h+ u PlKIL 2 (n) < -IteA h (v h ,v h ) + 2\v h \ lih + —Jo{v h ,v h ) 


7o 


1 2 " C 


< — Re A h (v h , v h ) + -u p|K||£ 2(n) + —11v-fc111,2 (an) 

C C 1 

H-Jo(v/,,v h ) + 9 —JiK,v h ). 

7o u-ph 2 7 i 

Thus, subtracting both sides of the above inequality by 5 W 2 p||v^||| 2fn) and using 
both (2.6) and (4.7) yield 


MU + w 2 p||v h ||| 2(n) 

< —2Re4(v /l ,v /l ) + 


1 


1 


<c(i + — + 


1 


+ 


7o uihcA uj 2 ph 2r yi 


1 

7o uihcA ' u} 2 ph 2r Yi 
'\A h {v h ,v h )\. 


Im A h (v h , v h ) 


Hence, (4.4) holds. The proof is complete. □ 

Remark 4.5. ( |4.4| ) -(4.5) is called weak coercivity because of the complex magni¬ 
tude (instead of the real part) is used in the left-hand sides of these inequalities. 

As a consequence of the discrete weak coercivity of Theorem |4.4| we readily obtain 
the following stability estimates. 

Theorem 4.6. Suppose that u h G V/, solves the IP-DG method given by (3.6). 
Then the following inequalities hold: 


(4.8) \u h \i th + UJ- p\\u h \\ 2 L 2 {n) < C'^C 2 ta ||f||| 2 (Q) + C'stallgll i,2( F ) y 

Cf \ 

(4.9) J 0 (u h ,u h ) + Ji(u h ,u h ) + (du ft ,u A ) r - ~(Csta||f||z,2(Q) + ||g|||2( r )J, 
where 


(4.10) 


C* s ta — — 


1 


1 


oj co 2 h ijj 3 h 2 71 ’ 

and C is a positive constant independent of u>, h, 70 , and 71 . 
Proof. By (4.4) and (4.51 we get 


£ — H—j 
7o 


KIU + W 2 


Pll u ^lll 2 ( n) + w ( 1 + ^ + ^ + UJ ip h ‘ 2 ll )( Auh ’ Uh ^ 


- C ( 1 + ^ + ^ + ^W) 


|A /l (u/ l ,u/ l )| 


< C(l + — + (Jh + u 2 pK 2 ^ x ) (ll f H i2 Wll u ^IU 2 (n) + l|g|U 2 (r)||u /l || i2(r) ^ 
C 


;( 


1 + 1 + A + 


- u 2 pV ' 7o ' coh ' u*pV^) l|f|li2 ( n ) ” l ^ _ ll u,i Hi 2 (n) 

C / i 1 1 v, 1|2 

+ a ^( 1 + ^o + ^ + ZfWlJ IISlli2(r) 

+ u,CA ( 1 + ^ 0 + fk + a Pph^n) l|u/lll ^ 2 ( r )‘ 
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Substituting (2.6) into the above inequality infers ( |4.8 l. 
Now, combining ( |4.5| ) with (4.8) yields 

■h 


bining ( |4.5| ) with (4.8) yields 
r o(u/i,u h ) + Ji(uft,,u ft ) +w(Au ft ,u /l ) r 

< \A h (u h , u h )| 

|g|!L 2 (r)||u/ 1 ||L2 (r) 


< | Ah (ll/j, Ufo) | 

< ||f||L 2 (n)||u/ l ||L2(Q) - 

< “ II f II r2CO'i ^ 11 f 11 Z.2 (f2) d Csta||g||l2(r) 

P C A / 

2 


1 

UJp 2 


WCa I. ||2 
■ 2 ll u ^llz, 2 (r) 


+ a^4 l|g|l ^ + 2 

< Csta||f||| 2 (n) d llg|li 2 (r) d — I u a|| z,2(r)- 

Up v ' UCA 2 v ’ 


Using (2.6) in the above inequality infers the desired estimate (4.9). □ 


Remark 4.7. When h is restricted to the pre-asymptotic mesh regime character¬ 
ized by the constraint u^h = 0(1) for 1 < /? < 1 + a, then we have 


(4.11) 


O s ta < C (-b U‘ 


c , ,2/3—3 

S I . ./3—2 , ^ 


7l 


Thus, the constant in the above stability estimate can be replaced by one that is inde¬ 
pendent of h (but still depends on Q) in the pre-asymptotic mesh regime. 

We conclude this section by establishing the following existence and uniqueness 
theorem. 

Theorem 4.8. For any u,h, 70,7i > 0, f £ L 2 (fl), and g £ L 2 (<9fl) there exists 
a unique solution u; t of (3.6). 


Proof. Due to the linearity of the PDE, it is easy to check that (3.6) is equivalent 


to a linear (algebraic) system. The stability estimate (4.8) immediately implies that 
the linear system only has the trivial solution when f = 0 and g = 0, which then 
infers that the coefficient matrix of the linear system is nonsingular. Thus, the system 


must have a unique solution and consequently (3.6) has a unique solution. □ 


5. Error estimates. The goal of this section is to provide error estimates that 
hold in all mesh regimes, in particular, in the pre-asymptotic mesh regime. In order 
to accomplish this we shall adapt a non-standard argument of sum by defining an 
elliptic projection operator onto and establishing estimates for the error between 


the solution u £ EDH d (U) of (1.1)-(1.2) and its elliptic projection, and utilizing the 
unconditional stability estimate to control the error between the projection and the 
IP-DG solution. In this section we assume u £ H 2 (fl) satisfying the estimate given 
in Theorem 14.11 

First, we show that the sesquilinear form a/ 1 (-, •) is continuous and satisfies a weak 
coercivity condition. 

Theorem 5.1. For any v,w £ E there exists a positive constant C independent 
of u, h, 70 , 71 such that 

( 5 . 1 ) K(v,w)| < cjlMiii^iiiwiiii^. 

Also for any 0 < S < 1 and £ Vj, there holds 

(5.2) Rea^v^v^) + ^1 - S + Im a h (v h ,v h ) > (1 - 5)||v^||f h . 
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Proof. (5.1) is easy to prove, we omit it and leave it to the interested reader. To 


show ( |5.2| ), we note that 

Rea h (vh,v h ) > K|^ + 2 Y |([v ft ], Mv h )n e }) e 


ee£l 


> Whllh - 2 ^ ||[v /l ]|| i 2 (e) ||{cr(v /t )n e } 


I L*{e) 


ee£l 


> 


\vh\lh- Y ^7|llM|| L 2 (e) (lkK)|| i2(ife) + || ( T(v^|| L2(if , ) ) 

ee£/ ™ "e 


> (l-£)|v h |f —Jo{vh,v h ). 
7o 


Also, 


Im a h {v h ,v h ) = J 0 (v h ,v h ) + Ji(v h ,v h ). 


(5.2) follows from combining the above two results. □ 


5.1. Elliptic projection and its error estimates. For any w £ E we define 
its elliptic projection wy £ Vj, as the solution to the following problem: 


(5.3) a h (w h ,v h ) + iu(Aw h ,v h ) r = a h (w,v h ) + iu(Aw,v h ) r Vv h eV t . 

Remark 5.2. By Theorem \5.1\ the sesquilinear form ah( •, •) + io;(A-,-)r is both 
continuous and coercive when 70 is taken to be sufficiently large. Thus, the Lax- 
Milgram Theorem ensures that wy is well-defined for such a choice ofjo. 

By the consistency of ah(-,-), it is trivial to check that the following Galerkin 
orthogonality property holds. 

Lemma 5.3. Suppose that w £ E and wy £ is its elliptic projection, then 


(5.4) a h {w - w h ,v h ) + iw(A(w - w h ), v h ) = 0 Vv h eV h . 


Now we are ready to estimate the error between the solution u of (1.1 H 1.2) and 
its elliptic projection ip,,. 

Theorem 5.4. Let u £ H 2 (fl) solve |Ll|)-([L2t and let £p £ V h be its elliptic 
projection defined in (5.3). Then there hold the following estimates: 


(5.5) 111u - UftIUpft, + w 2 £||u - u 7 l || i 2 (r) 

< Cf 2 h(£ + + ivh) 2 ^ ^||f||L 2 (n) + l|g||L 2 (r)) ■ 

(5.6) ||u — u^|| L 2 ( n ) < Cf"h 2 {f + 71 + Loh ) (||f||i 2 (n) + ||g||L 2 (r))> 


where £ = 1 + y ^ 1 and C is a positive constant independent of uj, h, 70 , and 71 . 

Proof. Let ip £ Vh denote the Pi-conforming finite element interpolant of u on 
%■ Then the following estimates are well-known (c.f. Hi): 

(5-7) ||u-u /l || L 2 (n) <Ch 2 \u \h 2 (q), 

(5-8) ||V(u- u,01^2(0) <Ch\u\ H 2 {Q) . 
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By the trace and inverse inequalities we get 

(5-9) 11|u — Uhllli,/* < C'(£ + 7i)5/i|u| ff 2( f2 ), 

(5-10) ||u — u/ l || i 2 (r ) < Chi | u | ff 2(Q). 

Set il> h := u^ — u/j. By Galerkin orthogonality along with the fact that %p h + u — = 

u — u/i we obtain 

ah(tp h , ip h ) + i *l>h) r =a h (u- u h , rj) h ) + iw(,4(u - u h ),ip h ) r . 

Next, it follows from ( |3.8[ ) and (5.2) with S = ! 2 (noting that C\ > |) that 

A 111 11 ll,/t < C^W^hWlM 

(l Ci\ 

< (7£ Rea h (ip h , ij) h ) + C£ Ima h (ip h ,ip h ) 


= C£ Re (a h {i/> h , i/> h ) + i u(Atp hl ip h 

+ C£ Q + Im(a h ('ip h ,%p h ) + iu(A-ip h ,-ip h ) r j 

- C& (j + (Mh: i’h ) r 
= Re (a h (u-u h ,ij} h ) + iu(A(u-u h ),ip h ) r ^ 

+ C£ Im (a h (u - u h ,-ip h ) + iu(A(u - u h ),^ h ) r ^j 

~ C& ^ (Aip h , ip h ) r 

< C'ClIlV’hllli./illlu- uHK.h + Cw^||i/j^|| L 2 (r) ||u - u|| L 2 (r) - Cu^\\ip h \\l 2(:r) 


+ C£ (lll^llluftlllu— u|||i, ft + w||'j/’Ji 2 (r) ||u —u 

< C'^ 2 HI^/Jlli,h|||u-u|||i i/t + 2 Ca;C 2 ||u- u ft ||| 2(r) - jw£ 2 ||VbJU 2 (r) 

< C^lllu- Uh\\\l, h + jlW^hlWlh ~ +2C'wC 2 || u -u/ 1 ||l 2(r) . 


Substituting (5.9) and (5.10) into the above estimate gives 

Ill^Jli.h +^ 2 Hh\\li(r) < C'^lllu- Ufcllli,/, +^ 2 ||u- Uh||! 2(r) ) 

< C(VC£ + 7i)^‘'l u lff2(Q) + w^ 2 /i 3 |u|^2(q)^ 

= CZ 4 h 2 (C + 7i + wh) |u|^ 2 (n) 

< C{; 4 h 2 (t; + 7 i + uh) ^ (l|f|| 2 2 (n) + Uglier))• 


W\fPh\\\i,h+^H\\iPh\\Li(r) 

< C( 2 h {£ + 7i + uh) 2 + ^2 ^ (||f||L 2 (n) + !|g|| 2 2 (r)) ■ 


Thus 
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Recall that u — u/i = u u/j if> h . By the triangle inequality we get 

|||u- U/ 1 ||| 1 ,/ 1 +o;5^||u - Uh||i 2 (r) 

< 11 |U — U h \\\ hh + W5£||ll- U h || L 2 (r) + Hl^hllll.fc + U^\\lp h \\ L 2 {T) 

< C £ 2 h ( t ; + 7i + uh) 2 ^ u a + — ^ (||f||L 2 (f2) + l!g|li, 2 (r)) • 


Therefore, (5.5) holds. 

To obtain (5.6), we appeal to the duality argument by considering the following 
auxiliary PDE problem: 


—div (cr(w)) = u — U/i in fl, 

er(w)n — iww = 0 on T. 

It can be shown that there exists a unique solution w £ H 2 (fl) such that 
(5.H) l|w|| ff 2 ( n) < C||u - u h || L 2 (n) 

Let Wft, £ V/,. denote the Pi-conforming finite element interpolant of w over 7h- Then 
by testing the above auxiliary problem with u — we get 

l|u- fl ftllz, 2 (n) = -(u- Ufc,,divcr(w)) n 

= a h ( u - u h , w) + iw(A(u - u^), w) r 
= a h (u - u h , w - w^)iw(A(u - u/*), w - w h ) r 

< c(j||u - Uft|||l,h|||w - Wh|||l,fc. + w||u - U h || L 2 (r) ||w - 

< c((^ + 7 i)’^ll w lli/ 2 (a)lll u - Uhllli./i +w/i^||w|| J? 2 (n) ||u- u/,||z,2 (r) ) 

< <?((£ + 7 i)^ll u - Uft||i2 ( n)|||u- u/jHIi.fc + wft 5 ||u- Uft||z,2(n)||u- Uh|U 2 (r))- 

Hence, 


|u- u/,||L2(n) < Ch(£ + 7i + uh ) 2 (|||u - u h |||i,h + w 2 £||u - u^|| L 2 (n) ^ 
< C£ 2 h 2 (£ + 7l + uh) + ~^J (|| f || L 2 (uj) + l|g|lL 2 (T))- 


Thus, (5.6) holds. □ 


5.2. Error estimates for the IP-DG method. The goal of this subsection is 
to obtain upper bounds for the error e^ := u — u^,, where u £ H 2 (fl) solve ( [IT] ) ( fl~2j ) 


and is its IP-DG approximation defined by (3.6). We shall accomplish this goal by 
using the stability estimate from Theorem |4.6 


and the elliptic projection studied in 


Subsection 5.1 Subtracting (3.6) from (3.5) immediately infers the following Galerkin 


orthogonality for e^: 


(5.12) 


a h (e h ,v h ) - u 2 p(e h ,v h ) n + iu(Ae h , v h ) r = 0 Vv h £ V,, . 


Let £ Vh be the elliptic projection of u as defined in (5.3). Next we define 
i/> := u u h and (f) h := — u/,. Thus, can be decomposed as = ip + 4> h ■ By 
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Galerkin orthogonality given in (5.4) and (5.12), we have the following identity: 


(5.13) a h (4> h ,v h ) -u 2 p(cp h ,v h ) n + iaj(A<f> h ,v h ) r 

= -a h {il !>, v h ) + u} 2 p(ip, v h ) n ~ i u(Aip, v h ) T 
= u 2 p(ip,v h ) n Vv,, e V,,. 


In other words we obtain that <p h £ solves (3.6) with f = pip and g = 0. This 


allows us to bound <p h by using the estimates from Theorem |4.6| and Theorem 5.4 
Specifically, we have the next lemma. 

Lemma 5.5. Let u £ H 2 (fl) solve |Tt|)-([L2|), u/i be its IP-DG approximation, 
and u/! its elliptic projection. Then <p h := u/j — satisfies 


(5-14) \\(p h \\ lih +ujp^\\(p h \\ L 2 {n) 

< Cf 2 u 2 h 2 C s + 7i + [fiJ 01 + ^ 2 ) (l|f||i 2 (n) + l!g||L 2 (r)) ■ 


Proof. As mentioned above, (5.13) implies that <p h solves (3.6) with f = to 2 pip 
and g = 0. Thus, an application of Theorem 4.6 yields 


\\4>h\\l,h + Up 2 \\(p h \\L 2 (n) < — C'sta||w 2 p'0||L2(Q) 

P 2 

< Cuj 2 p^ Csta || "01| L 2 (£7) ■ 


An application of Theorem 5.4 immediately infers (5.14). □ 

We are now ready to derive error estimates for our IP-DG method. The next 
theorem is a consequence of combining Theorem |5.4| and Lemma |5.5[ 

Theorem 5.6. Let u £ H 2 (f2) solve (1.1) (1.2) and be its IP-DG approxi¬ 
mation. Then u satisfies the following estimates: 


(5.15) 11u - u/jIIi >h + wp 2 11u - u ft || L 2 (n) 

< Cf 2 h(f + 71 + uh) (u a + (||f|! L 2 (n) + ||g|| L 2 (r)) 

+ Ct; 2 U}h 2 (l + wCsta) (f + 71 + + ^J2 ) l|f|U 2 (0)i 

(5.16) ||u — Uft,|| L 2 (Q) 

< Cf~h 2 (l + coGsta.) (f + 71 + uih ) (u> a + ^2^ ^||f||i 2 (n) + ||g||L 2 (3fi)) > 


where C sta 
h , 7oi 7i- 


and £ are defined in (4.10), and C is a positive constant independent of 


Proof. Recall that e/j = u — u/, = ip + (p h and the desired estimates for ip and <fi h 
have already been established in Theorem |5.4| and Lemma |5.5[ respectively. These 
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estimates are combined in the following steps to obtain (5.15): 


ll e /i||i,h + up 2 ll e ^l|L 2 (n) 

< IH^HIi./i + wp5||i/,|| i2(n) + \\(/) h \\i,h +wp3||0J i2(n) 

< Cf 2 h(f, + 71 + uti ) 5 (u a + 2 ^||f ||z, 2 (0) + ||g|U 2 (r)) 

+ Cf 2 uh 2 (f + 71 + ujh) + —'j (||f|| L 2 (n ) + ||g||z, 2 (r>^) 

+ C^ 2 Lo 2 h 2 C s ta,{£, + 7i + wft.) (ui a + —^ ^j|fi|z 2 (n) + ||g||L 2 (r)) 

< Cf 2 h(£ + 7i + w/i) (w“ + (||f||z 2 (n) + ||g|U 2 (n) 

+ Cf 2 ujh 2 (\ + wCgta) (£ + 7 i + ojh) (uj a + ) (j|f|U 2 (fi) + l|g!U 2 (r)^ ■ 


Similarly, (5.16) is obtained by combining Theorem 5.4 and Lemma 5.5 □ 


Remark 5.7. (a) The above estimates are sub-optimal in general due to the 

dependence of C s t a on the mesh size h. However, when h belongs to the pre-asymptotic 
mesh regime (i.e. when co^h = 0(1) for 1 < f3 < 1 + a), C s t a can be bounded by a 
constant which is independent of h as explained in Remark \4■ i\ Therefore, the above 
error estimates become optimal (in h) in the pre-asymptotic mesh regime. 

(b) Estimates that are optimal in h in the asymptotic mesh regime can be found 
in 


6. Numerical experiments. In this section, we present some numerical tests 
to demonstrate key features of the proposed IP-DG method. In all our tests we choose 
= (—0.5,0.5) 2 C M 2 (i.e. the unit square in M 2 centered at the origin), along with 
the material constants p = p = A = 1, and penalty constants 70 = 10 and 71 = 0.1. 
For the sake of testing the exact error, f and g are chosen so that the exact solution 
to the elastic Helmholtz problem is u = ^yy[e lwr — l,e~ lujr — 1] T , where r = ||x|| 
denotes the Euclidean length of x. This simple problem along with the subsequent 
numerical tests are chosen to mirror those in [5! for the scalar Helmholtz problem. 
Some sample plots are given in Figures [6.2| and [6.3| These plots demonstrate how well 
the proposed IP-DG method can capture the wave with large frequency when using 
a relatively coarse mesh. 

To partition the domain H, a uniform triangulation Th is used. For a positive 
integer n, define 7i/ n to be a triangulation of 2 n 2 congruent isosceles triangles with 
side lengths 1/n, 1/n, and \/2/n. Figure 


6.1 


shows the triangulation 7i/io- 


The numerical tests in this section intend to demonstrate the following: 

• absolute stability of the proposed IP-DG method, 

• error of the proposed IP-DG solution, 

• pollution effect on the error when ojh = 0(1), 

• absence of the pollution effect when u 3 h 2 = 0(1), 

• performance comparison between standard FE and the proposed IP-DG method 
on the test problem. 


6.1. Stability. In this subsection, the stability of both the proposed IP-DG 
method and the Pi-conforming finite element method will be discussed. Let u EEM 
denote the Pi-conforming finite element approximation of u. Recall that the proposed 
IP-DG approximation is unconditionally stable, i.e. it is stable for all u>,h, 7o,7i > 0. 
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Fig. 6.1. The triangulation T\/\q. 



Fig. 6.2. Plot of || Re (u7,)|R2(fp for u = 50 and h = 1/70. Both a top down view (left) and a 
side view (right) are shown. 


This has yet been established for the Pi-conforming finite element approximation. In 
fact, the stability of the Pi-conforming finite element approximation is only known to 
hold when h satisfies uj 2 h < C. 

plots both 11u/j.and ||u J ' £ ' M ||i^ for /i = 0.05,0.01 and w = 1, 2, • ■ ■ , 
200. We observe that Hu/jHi,/, decreases in a smooth fashion as w increases. This 
smooth behavior of is indicative of the absolute stability of the IP-DG ap¬ 

proximation. On the other hand, we observe oscillations in ||u^ BM ||i i i l that occur 
when we vary uj. This oscillation is indicative of the instability of the Pi-conforming 
finite element method when h is too large. 


Figure 


6.4 


6.2. Error. In this subsection, the optimal order of convergence for the proposed 
IP-DG method will be demonstrated. The pollution effect will also be examined. From 
Theorems 5.6 and the error analysis for the asymptotic mesh regime carried out in m 
(c.f. Remark |5.7| ) we expect the error in || • ||i^ to decrease at an optimal order in both 
the pre-asymptotic and asymptotic mesh regimes. In other words, ||u—u/j.||i,ft. = 0(h) 
is expected. Figure [6~5] is a log-log plot of the relative error ||u—u^Hi^/HuHi^ against 
the value 1/h for frequencies uj = 5,10,20,30. From this plot, it is observed that the 
relative error decreases at the same rate as h, thus displaying the optimal order of 
convergence in the relative semi-norm. Also displayed in Figure [63] is the error when 
uj varies according to the constraint ujh = 0.25. From this figure it is observed that 
the error increases as uj increases under this constraint. This is due to the pollution 
effect on the error for the elastic Helmholtz problem. 

The pollution effect for Helmholtz-type problems is characterized by the increase 
in error as uj is increased under the constraint uh = 0(1). This effect is intrinsic to 
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Fig. 6.3. Plot of || Re (u^) 11 z .2 (o) f or u = 100 an d E = 1/120. Both a top down view (left) and 
a side view (right) are shown. 




Fig. 6.4. Plots of \\ || 1 and \\u. EEM 


Helmholtz-type problems (c.f. 1T2J). It is well-known that the pollution effect can be 


the relative error is plotted against w as h is chosen to satisfy different constraints. 
Under the constraints wh = 1 and uih = 0.5, the pollution effect is present and the 
relative error increases as w is increased. On the other hand, when u> 3 h 2 = 1 is used 
to choose the the mesh size h, the pollution effect is eliminated. 


eliminated if h is chosen to fulfill the stronger constraint ui 3 h 2 = 0(1). In Figure 


6.6 


6.3. Comparison between the IP-DG method and the FE method. In 

this subsection, we compare the performance of the proposed IP-DG method and 
the Pi-conforming finite element method. As stated previously, the proposed IP-DG 
method is unconditionally stable while the Pi-conforming finite element method is 
only shown to be stable when h satisfies ut 2 h = 0(1). With this in mind, one can 
anticipate that in the case that the frequency oj is large, the IP-DG method should 
perform better. 

_ || Re(u,,)|| i 2 (Q) and || Re(u// BM )|| L 2 (n ) are plotted for tj = 100 
.20,1/200 on a cross-section over the line y = x. In addition, 


In Figures 


6.7 


and h = 1/50, 1/ 

|| Re(u)||z, 2 (f 2 ) is plotted to measure how well the respective approximations capture 
the true solution. In Figure |6.7[ it is observed that u/, already captures the phase 
of u with h = 1/50 while not fully capturing the large changes in magnitude. On 


FEM 


the other hand, for h = 1/50, has spurious oscillations. In this case, u ;i 

also fails to capture the changes in the magnitude of the wave. In Figure [GTS} we see 
that for h = 1/120, U/j captures the phase and changes in magnitu de o f the wave 
very well while u EEM still displays spurious oscillations. In Figure 


6.9 


we see for 

h = 1/200, both methods capture the wave well. However, the IP-DG method cap- 
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Fig. 6.5. Log-log plot of the relative error for the IP-DG approximation measured in the broken 
H 1 -seminorm for different values of u. 




Fig. 6.6. Relative error of the IP-DG approximations measured in the broken H 1 seminorm 
computed for different pairs of lj and h combinations. 


tures the wave slightly better. These examples demonstrate that the IP-DG method 
approximates high frequency waves better than the standard finite element method 
does, in particular, on a coarse mesh. This is of great importance when memory is 
limited or one wishes to employ a multi-level solver such as multigrid or multi-level 
Schwarz space/domain decomposition methods. 


I 

- 

J 1 

^ VVVVvvvs^. 


y = x 



Fig. 6.7. The left plot is of || Re(u^)||^ 2 (^) (solid red line) vs. || Re(u)|| L 2 (Q) (dashed blue line) 
for h = 1/50. The right plot is of || (solid red line) vs. || Re(u)||£ 2 /Q) (dashed 

blue line) for h = 1/50. 
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Fig. 6.8. The left plot is of || Re(uh)|| j, 2 (q) (solid red line) vs. || Re(u)|| £/ 2 (f 2 ) (dashed blue line) 
for h = 1/120. The right plot is of || Re(u^’ E ' M )|| I/ 2 (Q) (solid red line) vs. || Re(u)||_^ 2 (rj) (dashed 
blue line) for h = 1/120. 




Fig. 6.9. The left plot is of || Re(u^)||^ 2 (^) (solid red line) vs. || Re(u)|| L 2 (^) (dashed blue line) 
for h = 1/200. The right plot is of || B.e(u^ EM )\\ L 2 ^ (solid red line) vs. || Re(u)||^ 2 (Q) (dashed 
blue line) for h = 1/200. 
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